
*Create data set of CBSAs for inclusion in analyses
foreach yr in 11 12 13 14 {
	
		*Find all CBSAs with at least 1,000 admits - excluding hospitals in the top/bottom 1% of prices and not counting patients going to outside option hospitals
		use "$data\Master file`yr'.dta", clear
		*don't count pts who go to outside option
		drop if out_option == 1
		*don't count pts who go to outlier hospitals
		drop if ltp1 == 1 & gtp99 == 1
		gen obs = 1
		collapse (sum) obs_`yr' = obs, by(pt_cbsa prod) fast
		drop if pt_cbsa == "00000"
		
		save "$temp\admit_counts_`yr'.dta", replace		

}

use "$temp\admit_counts_11.dta", clear
merge 1:1 prod pt_cbsa using "$temp\admit_counts_12.dta", nogen
merge 1:1 prod pt_cbsa using "$temp\admit_counts_13.dta", nogen
merge 1:1 prod pt_cbsa using "$temp\admit_counts_14.dta", nogen
gen tot = obs_11+ obs_12+ obs_13+ obs_14
drop if tot == .
keep if tot>3999

contract pt_cbsa prod
preserve
keep if prod == "PPO"
save "$output\analysis_areas_ppo.dta", replace
restore
drop if prod == "PPO"
save "$output\analysis_areas_pos.dta", replace		

* * * * * Calculate patient level desc stats here * * * * * 
use "$data\Master file11.dta", clear
append using "$data\Master file12.dta"
append using "$data\Master file13.dta"
append using "$data\Master file14.dta"
drop if ltp1 == 1 & gtp99 == 1
drop if out_option == 1
merge m:1 pt_cbsa prod "$output\analysis_areas_ppo.dta"
keep if _m == 3
save "$temp\pt_data_ppo.dta", replace
gen pt_count = 1
bysort aha_id: gen uniq_hosp = _n == 1
egen  = total(uniq_hosp)
bysort aha_id: egen pt_per_hosp = total(pt_count)
egen avg_admits = mean(pt_per_hosp)
egen sd_admits = sd(pt_per_hosp)
egen avg_drg_wt = mean(drg_wt)
egen sd_drg_wt = sd(drg_wt)
gen pt_count = 1
collapse (sum) pt_count hosp_count, by(age_cat female hosp_count avg_admits sd_admits avg_drg_wt sd_drg_wt) fast
export excel using "$exp\Summary stats_`d'.xls", sheet("Pt Stats PPO") sheetreplace cell(A1) firstrow(variables)

use "$data\Master file11.dta", clear
append using "$data\Master file12.dta"
append using "$data\Master file13.dta"
append using "$data\Master file14.dta"
drop if ltp1 == 1 & gtp99 == 1
drop if out_option == 1
merge m:1 pt_cbsa prod "$output\analysis_areas_pos.dta"
keep if _m == 3
save "$temp\pt_data_pos.dta", replace
gen pt_count = 1
bysort hosp: gen uniq_hosp = _n == 1
bysort hosp: egen hosp_count = total(uniq_hosp)
bysort hosp: egen pt_per_hosp = total(pt_count)
egen avg_admits = mean(pt_per_hosp)
egen sd_admits = sd(pt_per_hosp)
egen avg_drg_wt = mean(drg_wt)
egen sd_drg_wt = sd(drg_wt)
gen pt_count = 1
collapse (sum) pt_count, by(age_cat female hosp_count avg_admits sd_admits avg_drg_wt sd_drg_wt) fast
export excel using "$exp\Summary stats_`d'.xls", sheet("Pt Stats POS") sheetreplace cell(A1) firstrow(variables)

use "$temp\pt_data_ppo.dta", clear
append using "$temp\pt_data_pos.dta", replace
gen pt_count = 1
bysort hosp: gen uniq_hosp = _n == 1
bysort hosp: egen hosp_count = total(uniq_hosp)
bysort hosp: egen pt_per_hosp = total(pt_count)
egen avg_admits = mean(pt_per_hosp)
egen sd_admits = sd(pt_per_hosp)
egen avg_drg_wt = mean(drg_wt)
egen sd_drg_wt = sd(drg_wt)
gen pt_count = 1
collapse (sum) pt_count, by(age_cat female hosp_count avg_admits sd_admits avg_drg_wt sd_drg_wt) fast
export excel using "$exp\Summary stats_`d'.xls", sheet("Pt Stats ALL") sheetreplace cell(A1) firstrow(variables)


* calcuate WTP for analysis areas.
foreach yr in 11 12 13 14 {
foreach prd in POS PPO {
	
	use "$output\analysis_areas_`prd'.dta", clear
	levelsof pt_cbsa, local(area)
	
	clear
	save "$temp\cbsa_`prd'_`yr'_data", replace emptyok
	
	clear
 	save "$temp\summary_`prd'_`yr'_data", replace emptyok

	*Calculate WTP measures for each CBSA identified above
	foreach a of local area {

		use "$data\Master file`yr'.dta", clear
		keep if prod == "`prd'" 
		drop if ltp1 == 1 & gtp99 == 1
		gen obs = 1
		
		*Flag CBSA for analysis based on patients CBSA - flag all observations from a given CBSA
		gen flag = pt_cbsa == "`a'"
		
		*Remove flag for obs associated with "outside option" hospitals
		replace flag = 0 if out_option == 1
		
		*Set CBSA flag to 0 for all obs associated with all hospitals with any admissions as the outside option
		bysort aha_hnpi: egen hosp = max(flag)
		replace hosp = 0 if out_option == 1

		*Flag all observations from any patient zip codes with admissions to hospitals with admits from the analysis CBSA
		bysort pt_zip: egen region = max(hosp)
		keep if region == 1
		
		*create a single "outside option" hospital for all hospitals that are too small/out of range/etc.
		replace aha_hnpi = "out" if out_option == 1 
		replace sys_hosp = "out" if out_option == 1 
		replace hosp_price = . if out_option == 1
		*This keeps all observations in the WTP calculation but the outside option WTP is not used in the regressions 
		
		save "$temp\reg input`yr'.dta", replace

		*Set parameters for groupings 
		local minsize = 25
		local chars = "pt_fips pt_zip mdc er_flag drg_type`yr' drg_wt_qnt drg age_cat female"
		local numchars = wordcount("`chars'")
		gen final_group = .
		gen group_size = .
		gen group_chars = .
			
		while `numchars' >= 1 {
			egen c = total(obs), by(`chars')
			egen g = group(`chars') if c >= `minsize' & final_group == .
			
			*define group number based on prior group num
			qui sum final_group
			local maxgroup = r(max)
			if `maxgroup' == . local maxgroup = 0
			replace final_group = g + `maxgroup' if final_group == .
			replace group_size = c if final_group != . & group_size == .
			replace group_chars = `numchars' if final_group != . & group_size !=. & group_chars ==. 
			
			*update char list 
			local tempchars = "`chars'"
			local numchars = `numchars'-1
			local chars = ""
			forvalues X = 1/`numchars' {
				local chars = "`chars' " + word("`tempchars'", `X')
			}
				
			drop c g
			replace obs = 0 if final_group != .
		}
			
		sum group_size, det
		sum obs if final_group == .
		drop if final_group == .
			
		*collapse groups
		replace obs = 1

		preserve
			collapse (sum) patients = obs, by(final_group group_chars group_size) fast
			gen analysis_area = "`a'"
			gen yr = `yr'
			gen prod = "`prd'"
			append using "$temp\summary_`prd'_`yr'_data"
			save "$temp\summary_`prd'_`yr'_data", replace
		restore

		collapse (sum) patients = obs, by(final_group aha_hnpi sys_hosp group_size) fast

		egen N_in_group = total(patients), by(final_group)
			
		*check group sizes
		assert N_in_group == group_size
		assert N_in_group >= `minsize'
		drop group_size

		*Patients = patients at the hospital in the group
		gen share = patients/N_in_group //share of total patients at the hospital 
		egen syspatients = total(patients), by(final_group sys_hosp) //total patients at the system
		gen sysshare = syspatients/N_in_group //share of patients at the system
			
		replace share = .95 if share==1
		replace sysshare = .95 if sysshare==1
			
		gen hospwtp = log(1/(1-share))
		gen syswtp = log(1/(1-sysshare))
			
		replace hospwtp = hospwtp * N_in_group
		replace syswtp = syswtp * N_in_group
			
		*calculate system WTP across all system hospitals
		preserve
			bysort final_group sys_hosp: keep if _n == 1 
			collapse (sum) syspatients syswtp, by(sys_hosp) fast
			save "$temp\sys_wtp_20`yr'.dta", replace
		restore
			
		*calculate hosp WTP & combine with sys WTPs
		collapse (sum) hospwtp patients, by(aha_hnpi sys_hosp) fast

		merge m:1 sys_hosp using "$temp\sys_wtp_20`yr'.dta"
		assert _m == 3
		drop _m
			
		merge 1:m aha_hnpi using "$temp\reg input`yr'.dta"
		keep if _m == 3
		drop _m
		
		gen pts = 1
		egen area_pts = total(pts)
		gen out_pts = aha_hnpi == "out"
		egen area_out_pts = total(out_pts)
		drop if aha_hnpi == "out"
		drop out_pts 
		
		ren (drg_wt`yr' cbsa_`prd'_hhi `prd'_hhi) (drg_wt cbsa_ins_hhi ins_hhi)
		
		gen pct_oop = (coins + copay) / calc_allwd
		
		gen analysis_area = "`a'"
		
		collapse (sum) pts (mean) pct_oop, by(analysis_area yr prod aha_hnpi sysid sys_hosp syswtp syspatients hospwtp patients ///
		hosp_state hosp_cbsa teaching forprofit bdtot ipdtot admtot mcrdc mcddc ftres hosp_price hosp_tech pct_unins pct_lt200pl ///
		cbsa_ins_hhi ins_hhi area_pts area_out_pts) fast
	
		
		* * create variables for regressions * *
		
		gen wtp_pp = syswtp/syspatients
		
		gen non_cbsa = length(hosp_cbsa)==2
		
		replace mcddc = 1 if mcddc == 0
		gen pct_mcd = mcddc/admtot
		qui sum pct_mcd, det
		replace pct_mcd = r(p99) if pct_mcd>r(p99)
		
		replace mcrdc = 1 if mcrdc == 0
		gen pct_mcr = mcrdc/admtot
		qui sum pct_mcr, det
		replace pct_mcr = r(p99) if pct_mcr>r(p99)
		
		*1 tech is .0067 = 1/147
		replace hosp_tech = .0067 if hosp_tech == 0

		gen ln_price = ln(hosp_price)
		gen ln_wtp = ln(wtp_pp)
		gen ln_beds = ln(bdtot)
		gen ln_mcr = ln(pct_mcr)
		gen ln_mcd = ln(pct_mcd)
		gen ln_hosp_tech = ln(hosp_tech)
		
		append using "$temp\cbsa_`prd'_`yr'_data"
		save "$temp\cbsa_`prd'_`yr'_data", replace
	}
}
}


exit
